#ifndef ARRAY_EXPR_FUNC_H
#define ARRAY_EXPR_FUNC_H

#include <limits>

//! for operator*()
template <unsigned int I, class return_T, class func_T, class... tuple_arguments_T>
struct expr_arg_tuple_evaluate
{
	inline static return_T unpack ( func_T&& func, 
									std::tuple<tuple_arguments_T...>&& tuple_args )
	{
		return expr_arg_tuple_evaluate<I-1, return_T, func_T, tuple_arguments_T...>::unpack( 
			std::forward<func_T>(func), 
			std::forward<std::tuple<tuple_arguments_T...>>(tuple_args), 
			std::forward<decltype(*std::get<I-1>(tuple_args))>(*std::get<I-1>(tuple_args)) );
	}    
		
	template <class ...function_arguments_T>
	inline static return_T unpack ( func_T&& func, 
									std::tuple<tuple_arguments_T...>&& tuple_args,
									function_arguments_T&&... function_args )
	{
		return expr_arg_tuple_evaluate<I-1, return_T, func_T, tuple_arguments_T...>::unpack( 
			std::forward<func_T>(func), 
			std::forward<std::tuple<tuple_arguments_T...>>(tuple_args), 
			std::forward<decltype(*std::get<I-1>(tuple_args))>(*std::get<I-1>(tuple_args)),
			std::forward<function_arguments_T>(function_args)... );
	}                       
};

template<class return_T, class func_T, class... tuple_arguments_T> 
struct expr_arg_tuple_evaluate<0, return_T, func_T, tuple_arguments_T...> 
{
	template <class... function_arguments_T>
	inline static return_T unpack (func_T&& func,
								   std::tuple<tuple_arguments_T...>&& tuple_args,
								   function_arguments_T&&... function_args )
	{
		return func( function_args... );
	}          
};

//! for operator[](int i)
template <unsigned int I, class return_T, class func_T, class... tuple_arguments_T>
struct expr_arg_tuple_subscript
{
	inline static return_T unpack ( int i, func_T&& func, 
									std::tuple<tuple_arguments_T...>&& tuple_args )
	{
		return expr_arg_tuple_subscript<I-1, return_T, func_T, tuple_arguments_T...>::unpack(
			i,
			std::forward<func_T>(func), 
			std::forward<std::tuple<tuple_arguments_T...>>(tuple_args), 
			std::forward<decltype(std::get<I-1>(tuple_args)[i])>( std::get<I-1>(tuple_args)[i] ) );
	}    
		
	template <class ...function_arguments_T>
	inline static return_T unpack ( int i, func_T&& func, 
									std::tuple<tuple_arguments_T...>&& tuple_args,
									function_arguments_T&&... function_args )
	{
		return expr_arg_tuple_subscript<I-1, return_T, func_T, tuple_arguments_T...>::unpack( 
			i,
			std::forward<func_T>(func), 
			std::forward<std::tuple<tuple_arguments_T...>>(tuple_args), 
			std::forward<decltype(std::get<I-1>(tuple_args)[i])>( std::get<I-1>(tuple_args)[i] ),
			std::forward<function_arguments_T>(function_args)... );
	}                       
};

template<class return_T, class func_T, class... tuple_arguments_T> 
struct expr_arg_tuple_subscript<0, return_T, func_T, tuple_arguments_T...> 
{
	template <class... function_arguments_T>
	inline static return_T unpack (int i, func_T&& func,
								   std::tuple<tuple_arguments_T...>&& tuple_args,
								   function_arguments_T&&... function_args )
	{
		return func( function_args... );
	}          
};

//! for fast_read(int i)
template <unsigned int I, class return_T, class func_T, class... tuple_arguments_T>
struct expr_arg_tuple_fast_read
{
	inline static return_T unpack ( int i, func_T&& func, 
									std::tuple<tuple_arguments_T...>&& tuple_args )
	{
		return expr_arg_tuple_fast_read<I-1, return_T, func_T, tuple_arguments_T...>::unpack(
			i,
			std::forward<func_T>(func), 
			std::forward<std::tuple<tuple_arguments_T...>>(tuple_args), 
			std::forward<decltype(std::get<I-1>(tuple_args).fast_read(i))>( std::get<I-1>(tuple_args).fast_read(i) ) );
	}    
		
	template <class ...function_arguments_T>
	inline static return_T unpack ( int i, func_T&& func, 
									std::tuple<tuple_arguments_T...>&& tuple_args,
									function_arguments_T&&... function_args )
	{
		return expr_arg_tuple_fast_read<I-1, return_T, func_T, tuple_arguments_T...>::unpack( 
			i,
			std::forward<func_T>(func), 
			std::forward<std::tuple<tuple_arguments_T...>>(tuple_args), 
			std::forward<decltype(std::get<I-1>(tuple_args).fast_read(i))>( std::get<I-1>(tuple_args).fast_read(i) ),
			std::forward<function_arguments_T>(function_args)... );
	}                       
};

template<class return_T, class func_T, class... tuple_arguments_T> 
struct expr_arg_tuple_fast_read<0, return_T, func_T, tuple_arguments_T...> 
{
	template <class... function_arguments_T>
	inline static return_T unpack (int i, func_T&& func,
								   std::tuple<tuple_arguments_T...>&& tuple_args,
								   function_arguments_T&&... function_args )
	{
		return func( function_args... );
	}          
};

//! because lambda function can not bind to generic type
struct tp_unpack_callback_push
{
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ expr.push( dim ); }
	int dim;
};
struct tp_unpack_callback_pop
{
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ expr.pop( dim ); }
	int dim;
};
struct tp_unpack_callback_load_stride
{
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ expr.load_stride( dim ); }
	int dim;
};
struct tp_unpack_callback_is_unit_stride
{
	tp_unpack_callback_is_unit_stride( int dim ) 
	: dim( dim ), result( true ) {}
	
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ result &= expr.is_unit_stride( dim ); }
	
	int dim;
	bool result;
};
struct tp_unpack_callback_suggest_stride
{
	tp_unpack_callback_suggest_stride( int dim ) 
	: dim( dim ), result( std::numeric_limits<int>::min() ) {}
	
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ if( int v = expr.suggest_stride( dim ) > result ) result = v; }
	
	int dim;
	int result;
};
struct tp_unpack_callback_is_stride
{
	tp_unpack_callback_is_stride( int dim, int stride ) 
	: dim( dim ), stride( stride ), result( true ) {}
	
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ result &= expr.is_stride( dim, stride ); }
	
	int dim, stride;
	bool result;
};
struct tp_unpack_callback_can_collapse
{
	tp_unpack_callback_can_collapse( int outer_dim, int inner_dim ) 
	: outer_dim( outer_dim ), inner_dim( inner_dim ), result( true ) {}
	
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) 
	{ result &= expr.can_collapse( outer_dim, inner_dim ); }
	
	int outer_dim, inner_dim;
	bool result;
};
struct tp_unpack_callback_advance
{
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) { expr.advance(); }
};
struct tp_unpack_callback_advance_n
{
	template <class expr_T>
	inline void operator()( int i, expr_T& expr ) { expr.advance( n ); }
	int n;
};

/////////////////////////////////////////////
//		function expression 
/////////////////////////////////////////////

template <class func_T, class... expr_T>
struct expression_func
{
	typedef decltype( (std::declval<func_T>())( *std::declval<expr_T>()... ) ) return_T;	
	typedef std::tuple<expr_T...> expr_arg_tuple_T;
	
	expression_func( func_T func, std::tuple<expr_T...> exprs  )
		: expr_tuple( exprs ), func( func ) {}
	
	expression_func( func_T func, expr_T... expr  ) 
		: expr_tuple( expr_arg_tuple_T(expr...) ), func( func ) {}

	return_T operator*() { return std::move( 
		expr_arg_tuple_evaluate<sizeof...(expr_T), return_T, func_T, expr_T...>::unpack(
		std::forward<func_T>(func), std::forward<expr_arg_tuple_T>(expr_tuple) ) ); }
		
	///////////	expression interface ///////////////////////////////////////
	void push( int dim ) 
	{ 	tuple_for_each( tp_unpack_callback_push{ dim }, std::forward<expr_arg_tuple_T>(expr_tuple) ); }

	void pop( int dim ) 
	{ 	tuple_for_each( tp_unpack_callback_pop{ dim }, std::forward<expr_arg_tuple_T>(expr_tuple) ); }
				
	void load_stride( int dim ) 
	{ 	tuple_for_each( tp_unpack_callback_load_stride{ dim }, std::forward<expr_arg_tuple_T>(expr_tuple) ); }

	bool is_unit_stride( int dim ) 
	{
		tp_unpack_callback_is_unit_stride cb( dim );
		tuple_for_each( cb, std::forward<expr_arg_tuple_T>(expr_tuple) );
		return cb.result;
	}
	
	int  suggest_stride( int dim ) 
	{ 
		tp_unpack_callback_suggest_stride cb( dim );
		tuple_for_each( cb, std::forward<expr_arg_tuple_T>(expr_tuple) );
		return cb.result;
	}
	
	bool is_stride( int dim, int stride ) 
	{
		tp_unpack_callback_is_stride cb( dim, stride );
		tuple_for_each( cb, std::forward<expr_arg_tuple_T>(expr_tuple) );
		return cb.result;
	}
	
	bool can_collapse( int outer_dim, int inner_dim ) 
	{ 
		tp_unpack_callback_can_collapse cb( outer_dim, inner_dim );
		tuple_for_each( cb, std::forward<expr_arg_tuple_T>(expr_tuple) );
		return cb.result;
	}
	
	void advance() 
	{	tuple_for_each( tp_unpack_callback_advance(), std::forward<expr_arg_tuple_T>(expr_tuple) ); }
	
	void advance( int n )
	{	tuple_for_each( tp_unpack_callback_advance_n{ n }, std::forward<expr_arg_tuple_T>(expr_tuple) ); }
	
	return_T operator[]( int i ) 
	{ 
		return std::move( 
		expr_arg_tuple_subscript<sizeof...(expr_T), return_T, func_T, expr_T...>::unpack(
		i, std::forward<func_T>(func), std::forward<expr_arg_tuple_T>(expr_tuple) ) ); 
	}
	
	return_T fast_read( int i ) 
	{ 
		return std::move( 
		expr_arg_tuple_fast_read<sizeof...(expr_T), return_T, func_T, expr_T...>::unpack(
		i, std::forward<func_T>(func), std::forward<expr_arg_tuple_T>(expr_tuple) ) ); 
	}
	
	constexpr typename bounding_checker<lbound_cb, expr_T...>::return_T 
	lbound( int dim )
	{ return bounding_checker<lbound_cb, expr_T...>::get_value( dim, expr_tuple ); }
	
	constexpr typename bounding_checker<shape_cb, expr_T...>::return_T
	shape( int dim )
	{ return bounding_checker<shape_cb, expr_T...>::get_value( dim, expr_tuple ); }
	
	constexpr typename bounding_checker<ordering_cb, expr_T...>::return_T
	ordering( int dim )
	{ return bounding_checker<ordering_cb, expr_T...>::get_value( dim, expr_tuple ); }
	////////////////////////////////////////////////////////////////////////
	
private:
	expr_arg_tuple_T expr_tuple;
	func_T func;
};

template <class func_T, class... Arg_T>
auto create_expression_func( const func_T& f, const Arg_T&... v ) ->
expression< expression_func<func_T, decltype(as_expression<Arg_T>::create(std::declval<Arg_T>()))...> >
{ 
	typedef expression_func<func_T, decltype(as_expression<Arg_T>::create(std::declval<Arg_T>()))...> expr_func_T;
	expr_func_T expr( f, as_expression<Arg_T>::create(v)... );
	return expression< expr_func_T >( expr );
}

template <class func_T>
struct expr_func_adapter
{
	expr_func_adapter( const func_T& func ) : func_(func) {}
	
	template <class... Arg_T>
	inline auto operator()( const Arg_T&... v ) -> 
	decltype( create_expression_func( std::declval<func_T>(), std::declval<Arg_T>()... ) )
	{ return std::move( create_expression_func( func_, v... ) ); }
	
private:
	const func_T& func_;
};

/**
	Example:
	auto func = []( double a, double x, double b ) { return a*x+b; };
	auto func_adp = create_expr_func_adapter(func);
	func_adp( any expressions... ); i.e.,
	func_adp( 2.0, array_type , 1.0 );
	func_adp( 2.0, array_type , array_type );
	func_adp( 2.0, array_type op array_type , 1.0 );
	...
 */
template <class func_T>
expr_func_adapter<func_T> create_expr_func_adapter( func_T&& func )
{
	return std::move( expr_func_adapter<func_T>( std::forward<func_T>(func) ) );
}

namespace spx
{
	#define DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( FUNC, FUNC_ADP )										\
																									\
	/** a wrapper to recover function from any expressable types to pure numeric function			\
		i.e., 																						\
		function_promoter< array<double, 1>, double > {												\
			double operator()( double, double )														\
		}																							\
		Arg_T can be any expressable types: array, expression, or primitive types					\
		function_recover provides operator () overloading to recover the function prototypes by 	\
		using expression evaluation operator *														\
	*/																								\
	template <class... Arg_T> 																		\
	struct function_recover_ ##	FUNC																\
	{																								\
		typedef decltype( 																			\
			FUNC_ADP ( *as_expression<Arg_T>::create(std::declval<Arg_T>())... ) 					\
		) result_type;																				\
																									\
		constexpr result_type operator()															\
		( decltype(*as_expression<Arg_T>::create(std::declval<Arg_T>()))... x )						\
		{																							\
			return FUNC_ADP ( x... );																\
		}																							\
	};																								\
																									\
	/* Arg_T can be any expressable type: array<T,DIM>, or expression<expr_T> */					\
	template <class... Arg_T>																		\
	auto FUNC ( const Arg_T&... v ) ->																\
	decltype( create_expression_func( std::declval<function_recover_ ##	FUNC<Arg_T...>>(), v... ) )	\
	{																								\
		return create_expression_func( function_recover_ ##	FUNC<Arg_T...>(), v... );				\
	}
	
	//! Trigonometric functions
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( cos, std::cos )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( sin, std::sin )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( tan, std::tan )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( acos, std::acos )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( asin, std::asin )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( atan, std::atan )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( atan2, std::atan2 )
	//! Hyperbolic functions
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( cosh, std::cosh )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( sinh, std::sinh )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( tanh, std::tanh )
	//! Exponential and logarithmic functions
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( exp, std::exp )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( ldexp, std::ldexp )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( log, std::log )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( log10, std::log10 )
	//! Power functions
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( pow, std::pow )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( sqrt, std::sqrt )
	//! Rounding, absolute value and remainder functions
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( ceil, std::ceil )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( fabs, std::fabs )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( floor, std::floor )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( fmod, std::fmod )
	//! Complex
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( abs, std::abs )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( arg, std::arg )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( norm, std::norm )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( conj, std::conj )
	DECLARE_GLOBAL_MATH_FUNC_AS_EXPR( polar, std::polar )
}

#endif //!< ARRAY_EXPR_FUNC_H